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Abstract  The  solution  of  certain  Toeplitz  linear  systems  is  considered  in  this  paper.  This  kind 

is  '~- 

of  system^a*©  encountered  when  we  solve  cert^  partial  differential  equations  by  finite  difference 
techniques  and  approximate  functions  using  higher  order  splines.  The  methods  presented  here  are 
more  efficient  than  the  Cholesky  decomposition  method  and  are  based  on  the  circulant  factorization 
of  the^*baaded  circulant"^matrix,  the  use  of  the  Woodbury  formula  and  algebraic  perturbation 
method.  ^ 
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1.  Introduction 

We  wish  to  consider  the  linear  system  of  the  form 


(1.1)  Az  =  f, 

where  the  coefficient  matrix  is  an  nth  order  symmetric  banded  matrix  of  Toeplitz  form 

/  ao  ai  ...  ap  \ 


or  cyclic  form 
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X  =  (xi,  X2, . , . , x„)^  is  the  unknown  n-vcctor,  and  /  is  the  given  right  hand  side. 

This  class  of  linear  systems  occures  in  solving  certain  kind  of  boundary  value  problems  by 
finite  difference  techniques,  solving  biharmonic  equation  by  Fourier  method,  and  in  higher  order 
spline  approximation [2,  3,  4,  6,  10]. 

System  (1.1)  with  coefficient  matrix  of  form  (1.2)  cam  be  solved  by  band  Cholesky  decompo- 
sition[7|  or  by  Toeplitz  factorization [6].  Although  the  operation  counts  of  the  two  methods  are 
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about  the  same  the  later  oi;e  requires  less  storages.  If  the  system  has  coefficient  matrix  of  form 

(1.3) ,  then  the  Cholesky  decomposition  is  expensive,  and  the  circulant  factorization  presented  here 
is  more  favorable  in  terms  of  not  only  arithmetic  operations  but  also  storage  requirements.  The 
methods  presented  in  this  paper  are  based  on  the  fact  that  under  certain  condition  the  matrix  in 

(1.3)  can  be  factored  into  two  simpler  circulant  matrices,  and  the  corresponding  circulant  system 
may  then  be  solved  by  using  the  Woodbury  formula[8].  Furthermore,  the  banded  Toeplitz  matrix 
may  be  treated  as  a  perturbation  of  circulant  matrix,  and  Toeplitz  systems  can  be  solved  by  the 
combination  of  the  circulant  factorization  and  the  use  of  algebraic  perturbation  method [9]. 

In  §2,  we  will  describe  the  method  for  factoring  a  symmetric  banded  circulant  matrix  into  two 
circulant  matrices,  and  then  use  the  factors  to  solve  the  band  circulant  system  in  §3.  The  methods 
for  solving  band  Toeplitz  systems  will  be  studied  in  §4,  and  finally,  some  numerical  results  will  be 
given  in  §5. 

2.  Factorization  of  banded  Circulant  Matrices 

To  factor  the  banded  circulant  matrix  given  by  (1.3)  we  consider  the  real  polynomial  with  the 
elements  of  the  matrix  as  its  coefficients 

(2.1)  4>{z)  =  ap^-l - -I - 1- 
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the  characteristic  function  of  matrix  A.  Assume,  without  loss  of  generality,  that  Op  =  1.  We  have 
the  following  theorem. 

Theorem  2.1.  If  matrix  A  is  strictly  diagonal  dominant,  i.e.  |oo|  >  2(|ai|  H - 1-  jopl),  then  there 

exists  a  real  polynomial  l{z)  —  0o  +  0iz  H - 1-  Ppz’*,  |/?o|  >  1,  I3p  —  1,  with  all  roots  outside  the 

unit  circle  such  that  the  characteristic  function  <^(z)  can  be  factored  as 

(2.2)  <l>(z)  =  ~l(z)  Hz-^). 

Proof.  We  show  at  first  that  the  polynomial  ^(z)  has  no  root  on  the  unit  circle.  If  there  exists  a 
number  zq  on  the  unit  circle  which  is  a  root  of  the  equation 


(2.3)  Hz)  =  0, 

then  Zq  =  e’*  for  some  real  $,0  <  $  <  2^.  Substituting  zq  into  (2.3)  we  have 

oo  =  ~  j^ui  - ^  +  e"'***^  j 

=  -  2  [ai  COS0  4 - 1-  Op cosp0] . 
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It  follows  that 


laol  <  2(|ail+  -  -t-lap|), 

which  is  a  contradiction  to  the  assumption  of  the  theorem. 

We  now  note  that 

4>(z)  = 

and  (2.3)  is  a  reciprocal  equation[l].  Thus  if  zq  is  a  root  of  (2.3),  then  so  is  Zq^.  It  follows  that 
<f>{z)  has  p  pairs  of  roots  z[^\  z^\  such  that 

-(*^1  —  ^  it  —  1  2  » 

Sj  -  (fc)>  K  - 

■^2 

and  are  outside  the  unit  circle. 

Let 

(2.4)  ((*)=n  (*-»!*’)■ 

We  now  prove  that  l(z)  is  a  real  polynomial.  If  all  the  roots  are  real,  then  p{z)  is  real;  if  some 
of  the  2i*^’s  are  complex,  then  their  conjugate  complex  numbers,  which  are  outside  the  unit  circle 
too,  are  the  roots  of  (2.3)  since  the  coefficients  of  the  equation  are  real.  So,  it  is  obvious  that  l{z) 
is  a  real  polynomial  and  satisfies  (2.2),  and  the  proof  is  completed. 


It  is  easy  to  verify  that  the  corresponding  circulant  matrix  Ac  can  be  factored  as 


(2.5) 


where 
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To  compute  the  factor  l(z),  we  solve  the  equation  (2.3).  When  p  =  2  it  is  well  known[l,  5] 
that  the  roots  of  equation  (2.3)  are  given  by 


1  j - 

Pi  =2  »?i  +  V ’ 


P2=2  [»?i  -  V’>i  "‘*J  ' 

P3  |»?2  +  \Jv2~*  1 


P4=2  »12-  V»/2-4l  , 


_  1 
”  2  . 


Cl  +  y  af  -  4ao  +  8  , 


where 


[  »?2  =  2  -  y/fl?  -  4oo  +  Sj  . 

Having  computed  the  roots  we  choose  the  two  roots  the  absolute  values  of  which  are  greater  than 
1  as  and  and  form  the  coefficients  of  the  factor  l{z)  via 


/So  = 


1/02  =  1. 

When  p  is  greater  than  2  we  have  to  use  some  numerical  method,  for  example  the  Newton- 
Raphson  method,  to  solve  equation  (2.3),  and  then  use  the  relations  between  the  roots  and  coeffi¬ 
cients  to  calculate  the  factor  l{z). 

S.  The  Solution  of  Band  Circulant  Systems 

In  this  section  we  will  use  the  circulant  factorization  described  in  the  previous  section  to 
develop  a  method  for  solving  the  band  circulant  system 

(3.1)  AcX^f 


as  well  as  computing  the  inverse  of  banded  circulant  matrices. 

It  is  evident  that  the  system  can  be  solved  by  solving  following  two  systems 


and 


(3.3) 
where 

(3.4) 

Let 


L  = 


and 


L^x  =  y, 

d  =  M- 

f  00 

0P 

V  0p  ■■■  00  J 
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0p 


Then  L  can  be  written 


(o) 


L=L+  ’’  i?(0^  /p), 


where  Ip  is  the  pth  order  identity  matrix  and  O  the  (n  — p)-by-p  zero  matrix.  Using  the  Woodbury 
formuia[8],  the  inverse  is  given  by 

1-1 


i2-*-h(0^  Ip)L 


"O 


[QT  Ip)L-\ 


and  the  solution  of  (3.2)  is 


=  i?-‘  +  (0^  (O^  Ip)L-^d, 


or 


(3.5)  y  =  h-Wg, 

where  h  =  (hi, ha,  •  •  •  ,hn)^,  W  and  g  are  the  solution  of  the  following  equations,  respectively 


Lh^d, 


(3.6) 
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LW  = 


(Ip\ 

\o  y’ 


Bg  =  r, 


(3.9)  B  =  i2-‘  +  (0^ 

To  compute  the  pth  order  matrix  B,  we  first  solve  the  equation  (3.7).  Since  IV  is  the  first  p 
columns  of  L~^,  and  L  is  a  lower  triangular  Toeplitz  matrix  and  so  is  its  inverse,  W  is  uniquely 
defined  by  the  first  column  of  L~*,  which  is  the  solution  of  the  equation 


(3.10) 


Lit;  =  (1,0,..., 0)^, 


and  can  be  computed  with  0{pn)  operations. 

Denote  by  Wi,W2, .  < . ,  t^n  the  components  of  vector  tr,  then  we  have 


(3.11) 


(3.12) 


W2  Ull 


.W„  W„-l  ...  W„.p+l 


(O^  Ip}W^ 


p+l  ^n—p  •  •  •  ^'n—2p+i 


W„  Wn-l 


V>n-p+l 


which  is  also  a  Toeplitz  matrix. 

The  matrix  R~^  is  an  upper  triangular  Toeplitz  matrix  and  can  be  calculated  with  0(p*) 
operations.  Thus  B  is  Toeplitz,  so  solving  (3.8)  will  cost  O(p^)  operations.  Having  computed  B 


'*•  **•  •*.  •  *  “’n  •*-  •* 
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and  solved  the  equations  (3.6),  (3.7)  and  (3.8),  the  auxiliary  vector  y  can  be  found,  and  we  can 
then  solve  equation  (3.3)  in  a  similar  way.  Since 


and 


the  solution  vector  x  is  given  by 


(3.13)  x  =  r-Vs, 

where  r  =  (ri,  rj, . . . , r„)^  is  the  solution  of  the  equation 


(3.14) 

and 


L^r  =  y. 


f  ton-p+l 

V^'n-p+2  ■ 

Wn 

^n— p 

W-’fi— p+1 

Wn~l 

(3.15) 

Wi 

XV2 

W^n-p+1 

Wi 

and  8  is  the  solution  of  the  equation 


XVl  J 


(3.16) 


B^s=  {Ti,r2,...,Tpy  . 


The  asymptotic  operation  counts  of  the  method  would  be  0(5pn)  excluding  the  amount  of  work  to 
calculate  the  factor  l{z).  In  most  usual  case,  p  =  1  or  2,  and  finding  l{z)  does  not  cost  much  work. 
The  algorithm  may  be  summarized  as  follows. 

Algorithm  BCS  (Banded  Circulant  Solver)  solves  banded  circulant  system  (3.1).  -Assume  that  the 
parameters  . . .,  are  precomputed. 

1.  Solve  equation  (3.6)  for  h  by  forward  substitution. 

2.  Solve  equation  (3.10)  and  form  W  via  (3.11). 


3.  Compute  R~^  by  backward  substitution,  and  form  matrix  B. 

4.  Solve  equation  (3.8)  for  g  using  a  Toeplitz  type  method. 

5.  Calculate  the  solution  vector  y  of  (3.2)  via  (3.5). 

6.  Solve  equation  (3.14)  for  r. 

7.  Form  V  via  (3.15). 

8.  Solve  (3.16)  for  s. 

9.  Compute  the  solution  vector  x  via  (3.13). 

endalgorithm 

Algorithm  BCS  can  be  modified  to  compute  the  inverse  of  banded  circulant  matrix.  Since  Ac 
is  a  symmetric  circulant  matrix  its  inverse  A~^  is  also  a  symmetric  circulant,  which  is  uniquely 
defined  by  its  first  column,  that  is  the  solution  of  the  equation 


(3.17) 


A,u  =  (1,0,...,0)^. 


The  algorithm  BCS  may  directly  be  employed  to  solve  equation  (3.17).  But  in  this  case  the 
first  two  steps  of  the  algorithm  ate  essentially  the  same,  so  we  obtain  the  following  algorithm  for 
inverting  banded  circulant  matrix  requiring  0(4pn)  operations  by  modifying  the  first  two  steps  of 
the  algorithm  BCS  and  computing  the  solution  of  equation 


(3.18) 


Ly  =  /9o(l,0,...,0)^, 


instead  of  equation  (3.2)  in  step  5  of  the  algorithm  BCS. 

Algorithm  BCI  (Banded  Circulant  Inverse)  computes  the  inverse  of  banded  circulant  matrix.  As¬ 
sume  that  the  parameters  /?o,  ■  ■  ■•,  0p  are  precomputed. 

1.  Solve  equation  (3.10)  and  form  W  via  (3.11). 

2.  Compute  h  =  0ow. 

3.  Compute  i?"*  by  backward  substitution,  and  form  matrix  B. 

4.  Solve  equation  (3.8)  for  g  using  a  Toeplitz  type  method. 

5.  Calculate  the  solution  vector  y  of  (3.18). 

6.  Solve  equation  (3.14)  for  r. 

7.  Form  V  via  (3.15). 

8.  Solve  (3.16)  for  s. 


9.  Compute  the  first  column  of  the  desired  inverse  via  (3.13)  and  form  it. 
endalgorithm 


4.  Band  Toeplitz  Systems 

The  band  Cholesky  decomposition  is  an  efficient  method  for  solving  general  band  symmetric 
systems[7|,  and  it  can  of  cource  be  used  to  solve  band  Toeplitz  system 

(4.1)  Atx=f. 


But  the  application  of  this  method  to  Toeplitz  systems  not  only  costs  a  lot  of  arithmetic  operations 
but  also  requires  a  great  amount  of  storages  since  it  does  not  take  the  advantage  of  the  structure  of 
Toeplitz  matrix.  Fischer  etc. [6]  proposed  the  Toeplitz  factorization  method  for  the  solution  of  band 
Toeplitz  systems,  which  has  some  advantages  both  in  terms  of  arithmetic  operations  and  storage 
requirements.  In  this  section  we  will  use  the  circulant  method  described  in  last  section  to  develop 
an  alternative  to  the  Toeplitz  factorization  for  solving  band  Toeplitz  system  (4.1). 

Banded  Toeplitz  matrix  At  may  be  considered  to  be  a  (2p)'rank  perturbation  of  the  banded 
circulant  matrix  Ae,  i-e. 


(4.2) 

where 


A,  =  Ac- 


O 


\U{0^  /p)-  ,  ]u^{ip  o^). 


u  = 


ai\ 


Qp  ) 


Substituting  (4.2)  into  (4.1)  we  have 

(4.3)  /^)x-  ^^^£7^(4  0^)1  =  /. 

If  matrix  At  is  strictly  diagonal  dominant,  then  the  corresponding  circulant  matrix  Ac  is 
likewise,  and  therefore  is  nonsingular,  and  from  (4.3)  we  have 

O 


(4.4) 


'  [oj 


Ip)x-A, 


-1 

'C  ' 


Let  =  (ii,....Tp)^,  x^^'>  =  (ip+i,...,j„-p)^,  and  =  (i„_p+i, . . .  ,i„)^,  and 
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"-"'■O’ 

which  are  the  n-by-p  submatrices  consisting  of  the  first  and  the  last  p  columns  of  matrix  * , 

respectively.  Then  equation  (4.4)  becomes 

(4.5)  a:  =  y  + 

which  shows  that  the  solution  to  equation  (4.1)  is  the  linear  combination  of  the  solution  of  the 
corresponding  circulant  system 

(4.6)  Acy  =  / 

and  the  first  p  and  the  last  p  columns  of  the  inverse  of  the  corresponding  circulant  matrix. 

The  solution  to  (4.6)  can  be  obtained  by  algorithm  BCS  in  0{hpn)  operations,  and  the  inverse 
of  Ac  can  be  calculated  in  0(4pn)  operations  by  using  algorithm  BCI.  The  inverse  A~^  is,  as  we 
pointed  out  above,  symmetric  circulant  and  defined  by  its  first  column,  the  elements  of  which  are 
denoted  by  «i,  U2,  ...,  «n  satisfying 

Un— I  =  i  =  0, 1, . . . ,  [(n  “  2)/2j , 


where  [aj  is  the  integer  floor  function  of  z.  We  then  have 
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To  compute  the  first  p  and  the  last  p  components  of  the  unknown  vector  x,  we  premultjply 
equation  (4.5)  by  {Ip  O^)  and  (O^  Ip),  respectively,  resulting  the  following  linear  system 

f  {Ip  -  MipU^) 

I  +  {Ip  -  M[pU) 

where  Mn  and  A/ip  are  the  pth  order  submatrices  of  A~^  at  the  northwest  and  northeast  corner, 
respectively,  i.e. 


M„  =  . 


Mi„  = 


Wn— p+l 


«n-2p+2 


Un-p+l 


and  y(®)  are  the  p- vectors  with  the  first  and  the  last  p  components  of  vector  y  as  their  elements, 

respectively. 

Forming  the  coefficients  of  equation  (4.9)  will  cost  0(2p*)  operations  and  (4.9)  can  be  solved  by 
Gaussian  elimination  with  0(8p®)  operations.  Having  calculated  y,  u,  and  the  subvector 
x(2)  can  be  obtained  via  (4.5)  with  0(2pn)  operations.  When  p  <  n,  the  asymptotic  operation 
counts  of  the  algorithm  would  be  0(1  Ipn)  excluding  the  amount  of  work  to  compute  the  factor 
l{z).  The  algorithm  thus  proceeds  as  follows. 

Algorithm  BTS  (Band  Toeplitz  Solver)  solves  band  Toeplitz  system  (4.1).  Assume  that  the  pa¬ 
rameters  A),  /5i,  .. .,  0p  are  precomputed. 


1.  Solve  for  y  equation  (4.6)  by  using  algorithm  BCS. 

2.  Compute  the  first  column  vector  u  of  A~^  using  algorithm  BCI. 

3.  Form  and  solve  equation  (4.9)  for  and 

4.  Compute  vector  via  (4.5),  which  along  with  and  is  the  solution, 
endalgorithm 

5.  Numerical  Experiments 

The  algorithms  described  in  this  paper  were  tried  on  the  APVAX  of  the  Department  of  Com¬ 
puter  Science,  Yale  University,  and  compared  with  Toeplitz  factorization  and  Cholesky  decompo¬ 
sition.  The  programs  were  written  and  timed  in  FORTRAN. 

To  obtain  some  insight  of  the  accuracy  of  the  algorithms,  we  generated  a  number  of  vectors 
randomly,  which  were  considered  to  be  the  “exact”  solutions  and  then  multiplied  them  by  the 
coefiicient  matrices  to  generate  the  corresponding  right  hand  sides.  The  equations  were  solved  by 
using  the  algorithm  BCS  and  BTS  as  well  as  the  Toeplitz  factorization  and  Cholesky  method.  In 
all  the  experiments  the  results  differ  from  the  “exact”  solutions  only  in  the  last  digit,  indicating 
that  the  algorithms  presented  in  this  paper  are  stable. 

In  our  all  tests  we  let  p  «  2  and  chose  several  matrices  satisfying  the  assumption  in  theorem 
2.1.  The  execution  time  of  algorithm  BTS  and  the  Toeplitz  factorization  are  almost  the  same. 
For  solving  circulant  systems  the  algorithm  BCS  is  about  twenty  times  faster  than  the  Cholesky 
method  in  our  tests,  and  saves  a  lot  of  storages. 
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